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ABSTRACT 


Neutron stars are a remarkable marriage of Einstein’s theory of general relativity with nuclear physics. Their interiors harbor 
extreme matter that cannot be probed in the laboratory. At such high densities and pressures, their cores may consist 
predominantly of exotic matter such as free quarks or hyperons. Gravitational wave observations from the Laser Interferometer 
Gravitational-wave Observatory (LIGO) and from other interferometers, and X-ray observations from the Neutron Star Interior 
Composition Explorer (NICER), are beginning to pierce through the veil. These observations provide information about neutron 
star cores, and therefore, about the physics that makes such objects possible. In this review, we discuss what we have learned 
about the physics of neutron stars from gravitational wave and X-ray observations. We focus on what has been observed with 
certainty and what should be observable in the near future, with an eye out for the physics that these new observations will 
teach us. 


Key points: 


e The processes at play inside neutron stars encode general relativity, quantum mechanics, particle physics and nuclear 
physics that cannot be replicated in the lab. 


e Gravitational wave observations of binary neutron star mergers are beginning to provide information about the equation 
of state of supranuclear matter through constraints on the tidal deformability of neutron stars. 


e The X-rays emitted by hot spots on the surface of certain pulsars are beginning to provide information about nuclear 
physics through constraints on the radius of neutron stars. 


e Future observations of gravitational waves and X-rays from LIGO, Virgo, KAGRA and NICER will provide unprece- 
dented insights into the physics of neutron stars. 


Website summary: The observation of gravitational waves emitted in the merger of neutron stars and the observations of 
X-rays emitted by hot spots on their surfaces are beginning to reveal nuclear physics insights about these compact objects. 


1 The marvelous neutron star 


Neutron stars are wonderfully marvelous. With masses comparable to our Sun’s, but radii of only approximately 12 km, they 
are some of the most compact, and thus gravitationally powerful objects in the universe!. With monstrous magnetic fields”, 
sometimes a trillion times stronger than that of a refrigerator magnet, they funnel photons into beams that travel astronomical 
distances. With astounding rotation speeds that can reach up to hundreds of Hertz, rivaling professional kitchen blenders, they 
whip these magnetic fields and beams around, creating astrophysical lighthouses. Every time the beams cross Earth, radio 
telescopes record a pulse, and the counting of these pulses creates a clock. Neutron stars are, in fact, one of the most stable and 
accurate clocks known in the natural universe because of their marvelously stable rotation rates?. 

Neutron stars are also wonderfully unavoidable. A massive star is supported against gravity during most of its life by the 
radiation force that it produces in the thermonuclear fusion of its component gases. The by-products of this reaction are ever 


heavier elements, and thus, at some point in the star’s lifetime, it is no longer energetically favorable to continue the burning 
sequence. The massive star then sheds its outer layers in a supernova explosion, leaving behind an iron core that can no longer 
burn, and thus can no longer prevent gravitational collapse. In the iron core goes, just as dictated by Einstein’s theory of general 
relativity, and so if nothing came to the rescue, one would expect the formation of a black hole. 

Neutron stars, however, are wonderfully quantum mechanical. As the iron core contracts, it becomes energetically favorable 
for electrons and protons to combine to form neutrons and emit neutrinos via inverse beta decay. The iron core has now become 
essentially a soup of neutrons, peppered of course with a few other particles. But neutrons are fermions, and this soup is so 
dense that the distance between fermions becomes really small, forcing each neutron to feel the influence of its neighboring 
neutrons. By the Pauli exclusion principle, this then leads to a very large (Fermi) momentum and energy, and thus, a very large 
“degeneracy” pressure, whose gradient halts the collapse. 

As a result, neutron stars are wonderfully dense. Although their atmospheres (which comprise only the few centimeters 
closest to the surface) can have atoms (albeit possibly distorted into near-cylinders by the strong magnetic field), just a bit 
deeper, matter is so dense that electrons do not belong to individual nuclei. At densities greater than ~ 10’ g cm~, the Fermi 
energy of electrons becomes high enough that the matter becomes progressively richer in neutrons*, which produces nuclei 
such as '?°Rb, which have 40 protons and 80 neutrons. At densities greater than 4 x 10!! g cm7? it becomes possible for 
neutrons to “drip out” of the nucleus, which means that matter is a mix of free neutrons, free electrons, and nuclei. At even 
higher densities, nuclei can cluster together to form “pasta-like” structures, such as one-dimensional strings (“spaghetti") or 
two-dimensional surfaces (“lasagna"). Going in a bit deeper, at about “nuclear saturation density” (2.7 x 10!4 g cm7?, so called 
because it corresponds to the density at the center of large nuclei), there are no longer any isolated nuclei, and we now have the 
neutron soup mentioned above, as shown in Fig. 1. Pushing to yet higher densities, it may be energetically favorable to form 
other baryons with at least one strange quark, such as hyperons ®6, until eventually close to the center of the star, quarks may 
become deconfined’ and one may encounter a degenerate quark-gluon plasma. 
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Figure 1. Cartoon structure of a neutron star and its internal structure. 


Neutron stars are therefore wonderfully cool, but not just because they are physically interesting: they are also “cold.” 
Imagine a neutron star composed of an immense number of fermions, which due to the Pauli exclusion principle must occupy 
different energy states. The highest energy state, the Fermi energy, is therefore enormous, and one can think of these fermions 
as forming a kind of “gas.” If so, you can associate an effective temperature to this gas, by dividing the Fermi energy by the 
Boltzmann constant. This temperature turns out to be remarkably high, on the order of 10! K, about 5 orders of magnitude 
higher than the temperature at the center of the Sun. Therefore, even though isolated neutron stars are actually very hot in 
terms of their actual temperature, which at their cores might typically be 108 — 10!° K, they are cold relative to their Fermi 
temperature. Of course, when neutron stars are born, they can be significantly hotter; for example, proto-neutron stars can have 
core temperatures of 10!! K. But such stars cool down very rapidly via neutrino emission, dropping by two orders of magnitude 
in just a thousand years, which is an extremely short time scale by astronomical standards. Thus, unless one is observing their 
birth, neutron stars are, for most purposes, cold objects relative to their Fermi temperature. 

Although it may not seem so from the above description of their complicated structure, there are some ways in which 
neutron stars are wonderfully simple. Many macroscopic aspects of neutron stars that have been or could be observed using 
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electromagnetic radiation or gravitational waves, including their mass, radius, and tidal deformability, depend only on the 
nature of gravity (here we assume general relativity) and the equation of state’, which determines the pressure given other 
quantities such as the energy density, temperature, and composition. As argued above, the temperature is too low to make a 
difference in isolated and “old” neutron stars, although temperature cannot be neglected in the merger of neutron stars. The 
composition is usually assumed to be the equilibrium composition, because there is enough free energy to transition to the 
ground state (unlike, say, at terrestrial densities, where >°Fe is the ground state of matter but there is insufficient energy to cause 
fusion to guarantee that matter reaches that state). Other complications are also thought to be ignorable in the description of the 
equation of state; for example, shear and bulk viscosity are believed to contribute negligibly to the equation of state (although 
they may be relevant in the dynamics of the merger of neutron stars’). Thus, to an excellent approximation, in neutron star 
cores the pressure depends only on the energy density, which is to say that the equation of state is barotropic. 


What is the equation of state, and how can it be inferred from observations of neutron stars? 


As indicated in the text, the equation of state represents a mapping between the pressure and other quantities, such as 
the energy density, temperature, and composition. But in neutron stars, it is typically thought that in the core the pressure 
p depends only on the energy density £, so p = p(€). To see how observations of neutron stars can help constrain the 
equation of state, note that the equation of hydrostatic equilibrium, for nonrotating (and thus spherical) fluid stars in 
general relativity, is the Tolman-Oppenheimer-Volkoff (TOV) equation! 
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Here r is the coordinate distance from the center of the star and m is the gravitational mass inside this radius, while G 
is Newton’s gravitational constant and c is the speed of light. The quantity € includes the rest-mass density, so in the 
Newtonian limit where p < £, £ is dominated by rest mass-energy pc’. In this same limit, 2Gm/c* <r, so the above 
equation reduces to the Newtonian hydrostatic equilibrium equation dp/dr = —p(Gm/r). Given an equation of state 
p(€) and a central density £e, one can integrate this equation, combined with the continuity equation dm/dr = 4ar7e/c*, 
to find the (total) mass M and radius R of a star. In the Newtonian limit, larger €< guarantees larger M, but in general 
relativity there is a central density £c max beyond which further increase leads to lower M. This corresponds to an instability, 
and means that in general relativity a given equation of state has a maximum stable mass. Thus, any measurement or 
constraint on the mass and radius of a star (or other quantities such as the mass and tidal deformability of a star, or the 
mass and moment of inertia of a star) can be compared with the predictions of a set of equations of state. 

In these neutron star structure equations, however, there is no information about the composition of the core. Any 
composition that yields a p(€) that is consistent with observations will do. Thus, although observations of neutron star 
mass, radius, tidal deformability, moment of inertia, etc. provide valuable constraints, they cannot by themselves tell 
us whether the core is mainly neutrons, or hyperons, or free quarks, or something else. Some additional information 
can in principle be obtained by measurements of the temperatures of neutron stars of a given age and mass (because 
temperatures depend on transport properties that in turn have some dependence on composition) or about the gravitational 
waves emitted in the (hot) merger of neutron star binaries. However, at this time, temperature and merger information is a 
bit too uncertain or unavailable to provide strong and reliable constraints. 


2 The equation of state puzzle 


The problem is then “simple”: solve for the equation of state using many-body quantum mechanics and quantum chromody- 
namics, and then use this in conjunction with the Einstein equations to predict the observable properties of the equation of state. 
Unfortunately, this is easier said than done. The Einstein equations part is not the problem. In fact, it is relatively straightforward 
to compute the observable properties of neutron stars in general relativity, once you are given an equation of state. The problem 
is quantum mechanical, and relates to the “sign problem"!® present in quantum chromodynamics calculations at nonzero baryon 
chemical potential. Although this theory has straightforwardly calculable predictions!! for matter at finite temperature and zero 
net baryon density (i.e., matter with almost equal numbers of baryons and antibaryons), current numerical methods become 
exponentially more costly when the net baryon density is large. As a result, it is not possible to compute the properties of 
neutron star core matter using current first principles approaches. 

For this reason, an abundance of models for the equation of state of matter at supranuclear densities and effectively zero 
temperature have been put forth'*. One approach is to solve quantum-chromodynamic-motivated models using the best 
microphysics possible given the constraints of the sign problem, see e.g.!*-!°. Another is to assume that we know the equation 
of state up to some threshold density (which is typically around nuclear saturation density) and then extrapolate to higher 
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densities using some (often parameterized) one-dimensional function, such as a piecewise polynomial!’. In this latter approach, 
the aim is to use nuclear physics experiments and astrophysical observations to constrain this phenomenological function, and 


then to study what the constraints imply for nuclear physics!®. 


The phenomenological nature of the second approach does not imply that interesting nuclear and particle physics is 
unimportant. On the contrary, the interactions between nuclei, neutrons and quarks can greatly influence the functional form 
of the equation of state, and therefore, both approaches attempt to include as much physics as possible. For example, if the 
matter in the core of neutron stars transitions into deconfined quarks, then under certain circumstances it is possible for the 
speed of sound of the fluid (the square root of the derivative of the pressure with respect to the energy density) to become very 
small or zero!?. This is called a first-order phase transition in the quantum chromodynamic phase diagram, which for old and 
isolated neutron stars is two dimensional (pressure versus density or chemical potential only). The speed of sound also cannot 
exceed the speed of light (the so-called “causal limit”), and at extremely high energy densities it is expected to approach square 
root of one third the speed of light (the so-called “conformal limit”)*°. The latter arises because at sufficiently high densities, 
the particles’ energy is dominated by their Fermi momentum, so they can be treated effectively as a relativistic gas. It is not 
currently known if the conformal limit applies at the densities expected inside the cores of neutron stars, or whether other phase 
transitions may be present. The equation of state puzzle then may be solved through observations, since certain quantities 
that are (directly or indirectly) observable are determined by the equation of state (see the box for more details). Equations of 
state are sometimes called “stiff” or “hard" when the slope of the pressure-energy density curve is large. On the other hand, 
a “soft" equation of state, for which the pressure increases slowly as the energy density increases, has a smaller maximum 
mass. Equations of state that contain first-order phase transitions around nuclear saturation density allow for stars with similar 
masses but different radii, which have been dubbed mass twins!. Figure 2 shows a schematic representation of stiff and soft 
equations of state, and equations of state with first order phase transitions, together with their respective mass-radius curves. 
Other interesting properties of nuclear/quark matter include heat capacity, superfluidity, viscosity, and shear modulus that may 
be probed through cooling’, glitches**°, and r-mode oscillations”®~’ of neutron stars. Other basic physics questions that 
one could address through neutron star observations include the role of three-body forces and the threshold for the appearance 
of hyperons (the “hyperon puzzle”). 


stiff 


mass twin 


first-order 
phase 
transition 
1 3 5 10 12 14 
nB / Nsat R [km] 


Figure 2. Schematic representation of equations of state showing the square of the sound speed against baryon number 
density normalized to nuclear saturation density (left) and their corresponding mass-radius curves for neutron stars (right). 
Generally, stiff equations of state have larger sound speed, maximum mass and radius than soft ones. Equations of state with 
first-order phase transition (energy densities at which g = 0) can produce mass twins if they are large enough and if they occur 
at densities around what is shown in the figure. 
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3 A new dawn for nuclear astrophysics 


Up until the 2010s, most astrophysical information about neutron stars came from observations of the light they emit. Perhaps 
the best known of these are the radio pulses emitted by neutron stars in a binary orbit with another compact companion, such as 
another neutron star or a white dwarf. The timing of these pulses allow the careful reconstruction of the binary orbits, including 
relativistic effects, which earned Hulse and Taylor?® a Nobel Prize in Physics in 1993. Among the many wonders these binary 
pulsars have revealed, what concerns us the most here is the inference of the masses of the binary companion. To date, the 
heaviest neutron star with a well-measured mass is PSR JO740+6620 (where the numbers represent the location of the source in 
the sky in astronomical notation) at m = 2.08 +0.07 Mo”. Since high neutron star masses can be compared with the maximum 
mass predicted using different equations of state, the observation of PSR JO740”? and similarly high-mass neutron stars*”*! 
rules out equations of state that are too soft. Indeed, neutron star mass measurements based on radio timing have been, and 
continue to be, a lynchpin for astronomical constraints on matter beyond nuclear density, because the data, analysis, and model 
inferences are all well-understood. 


A new window to the universe and to neutron stars opened in 2015 when the advanced Laser Interferometer Gravitational- 
wave Observatory (or advanced LIGO for short) made the first direct detection of gravitational waves**. These waves are 
perturbations to gravity predicted by Einstein and produced when massive objects accelerate. Because of the feebleness of 
the gravitational force, a tremendous amount of mass has to accelerate to a tremendous degree for us to be able to observe 
gravitational waves that originate a cosmological distance away from Earth. But this is exactly what happened on 14 September 
2015. The GW150914 event (named after its discovery date) was shown to be produced by the merger of two black holes with 
masses roughly 30 times that of our Sun, at velocities close to half the speed of light**. This single event announced the birth of 
gravitational wave astrophysics. 

Many discoveries have been made in gravitational waves since 2015*°, but the most exciting for us was made in 2017°4. On 
August 17th of that year, advanced LIGO (and this time, also another gravitational wave detector called Virgo, too) detected 
gravitational waves again, but this time the frequency at which the signal peaked was much higher than in 2015. For objects 
such as black holes and neutron stars, whose radius R is just a few times Gm/c?, a higher frequency is a tell-tale sign of a 
much lower-mass merger. This is because Kepler’s Third law tells us that the square of the orbital frequency, which is directly 
proportional to the square of the gravitational wave frequency, scales linearly with the total mass of the binary and inversely with 
the separation d cubed. Thus at merger, when d ~ R ~ Gm/c’, the gravitational wave frequency is inversely proportional to the 
binary’s total mass. A detailed analysis of the GW170817 event later revealed that it was likely produced by the coalescence of 
two neutron stars, with masses of ~ 1.3 — 1.4 Mz at a mere 40 Mpc (i.e., 130 million light years) away from Earth**. Ona 
human scale, that is a ridiculous distance (astronomical in fact!), but to an astrophysicist, this is incredibly close (it corresponds 
to a cosmological redshift of about 0.009.) The closer the event, the stronger the signal, so using this event it was possible for 
the first time to extract information about the equation of state from gravitational waves. 


533 


But how do gravitational waves carry information about the equation of state? When two neutron stars spiral into each other 
and collide, before the collision takes place they are tidally perturbed by each other’s gravitational field. Just like Earth acquires 
a tidal bulge due to the Moon, inducing high and low tides in the ocean, when a neutron star gets close to another compact 
object, it will tidally deform. This tidal deformation requires energy, and so the neutron star “borrows” it from the orbital 
energy, therefore forcing the binary to spiral in faster than it would have otherwise. A speed up in the rate of inspiral directly 
affects how the gravitational wave frequency changes with time, because, as we said before, the orbital and gravitational wave 
frequencies are linearly related. Therefore, by carefully monitoring the evolution of the gravitational wave phase, one can 
in principle extract information about how much the objects that produced the wave were tidally deformed on their way to 
coalescence? 36 (see Fig. 3). 

This is exactly what the LIGO Scientific and Virgo collaborations measured from the GW170817 event?438-3?, The signal 
to noise ratio was large enough that, from the gravitational wave data alone, a double neutron star merger in which both stars 
had the same equation of state was somewhat preferred over other options, including two black holes (which is also strongly 
disfavored by the subsequent electromagnetic emission) and one neutron star and one black hole. If the event involved a binary 
neutron star system similar in its spins to those we observe in the Milky Way, then the probability distribution of the tidal 
deformability for both stars peaks at a nonzero value and implies a radius between ~ 10.5 km and ~ 13.5 km for both stars, at 
90% credibility*®. 


4 From gravitational waves to dense matter 


But how on Earth do you constrain the radius of neutron stars from a measurement of the tidal effect on the gravitational waves 
emitted? In very broad terms, there are essentially two approaches that have been pursued to do so. 

The standard approach is to adopt some model for the equation of state!”4°~? and then carry out Bayesian parameter 
estimation. In general, all equation of state models contain some low-density part, which is then extrapolated in some specific 
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Figure 3. Gravitational waveform from the last few cycles of a compact binary inspiral with (cyan solid) and without (black 
dashed) the effect of tides. The former (latter) corresponds to waves generated by inspiraling black hole (neutron star) binaries 


with vanishing (non-vanishing) tidal deformabilities. The waveform data is taken from’. 


mathematical way to higher densities beyond some threshold density (which has been typically taken to be between half and 
twice nuclear saturation density). One then employs standard Bayesian inference: for a given draw of the equation of state 
model and two draws of the central densities, one predicts the masses of the binary companions and their tidal deformabilities, 
and from this, the gravitational wave model that one compares to the gravitational wave data. Based on the relative likelihoods 
of the various draws in the exploration of the likelihood surface, one can then obtain posterior distributions for the pressure 
versus energy density, the mass versus radius, or other quantities. 

The main caveat of this approach stems from the need to prescribe an equation of state model. Any given model will place 
greater prior weight on some portions of equation of state parameter space than on others. For example, some models may 
exclude the possibility of first-order phase transitions, whereas others might strongly emphasize them. Some models may 
exclude wiggles, kinks and other crossover-type structure in the speed of sound, while others might focus on them. A reasonable 
approach then is to use a few distinct models in the hope (which fortunately is more and more the reality) that the data are 
informative enough that distinct but reasonable models lead to similar posteriors on the mass, radius and other observables. 

Another complementary approach relies on relations (often called “universal relations") between neutron star properties that 
are insensitive to the equation of state**. Particularly tight relations of this type include the so-called I-Love-Q relations**#° 
between the moment of inertia (I), the tidal deformability or Love number (Love) and the quadrupole moment (Q). For the 
purpose of inferring the radius of neutron stars from gravitational wave data, the most important universal relation is the binary 
Love one*®*’ (see*® for a similar relation), which relates the anti-symmetric combination of the tidal deformabilities to their 
symmetric combination and the mass ratio of the binary, and the Love-C one**;*?, which relates the tidal deformability of a 
star to its compactness. The binary Love relations can be used to effectively reduce the number of independent parameters 
needed to describe the gravitational waveform, and thus, improve the precision with which the tidal deformabilities can be 
measured*® 47:9, With a posterior on the deformabilities, the Love-C relation then provides the compactness of each star, 
which when combined with the gravitational wave measurement of the binary component masses, yields the radius of each star. 


How do you use relations insensitive to the equation of state to learn about the latter? 
And what is the Love number anyways? 


Let us begin with the second question first. The Love numbers are a set of real numbers introduced by Augustus 
Edward Hough Love in the 1900s to describe the Earth tides caused by the Moon. But of course one can study the tides of 
any massive body caused by any external perturbation, including the tides of a neutron star due to its binary companion. 
The Love numbers are proportional to the tidal deformability, which describes how much a star deforms in response to an 
external perturbation?'~>*. More precisely, an external perturbation generically induces a redistribution of the isodensity 
contours inside a massive body, which in turn can be described through mass and current multipole moments of the mass 
distribution. The tidal deformabilities are then formally defined as the constants of proportionality that relate how much 
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of a multipolar deformations is induced in a star .4 due to an external tidal perturbation Yy, i.e. Mly = Ap P4. 

The external perturbation can be of two classes (even or “electric” and odd or “magnetic” parity), depending on 
how it transforms under parity. Each class, in turn, can be decomposed into multipole moments, with the leading-order 
perturbation produced by the electric-type quadrupole tidal tensor ¢;;. Given this, the tidal deformability can also be 
classified and decomposed in an analogous manner, with the leading-order tidal deformability being the electric-type 
quadrupole tidal deformability Ag po or just A for short. This deformability is then defined via the relation Q;; = 16j;, 
where Q;; is the induced mass quadrupole moment. The calculation of any tidal deformability requires the solution to the 
perturbed Einstein equations for a star that is being deformed by some “external universe,” such as a binary companion 
sufficiently far away. 

The tidal deformabilities of the neutron stars in a binary affect the gravitational waves they emit during their inspiral 
because they modify the orbital energy and the rate of gravitational wave emission. For an equal-mass binary, the 
tidal effect on the orbital energy changes the phase by 5.5 times as much as the tidal modification of the gravitational 
wave emission rate to the leading post-Newtonian order, irrespective of the equation of state. As described in the 
text, the tidal deformations modify the Hamiltonian of the binary system by adding a term of the form*>** 6H œ 
Ui2(A /m*)(my /m)>U}, +1 — 2, where U12 = Gm/r12 is the binary’s Newtonian potential, with m the total mass and 
rı2 the orbital separation, and A; is the electric-type, £ = 2 tidal deformability of star 1, which scales with its radius 
to the fifth power. Since a binary system composed of tidally deformed stars has a larger (i.e., less negative) binary 
binding energy, it takes less time for gravitational waves to drain this energy away, and for the binary to inspiral. Such a 
modification in the inspiral orbital dynamics imprints directly on the gravitational waves emitted. 

Gravitational wave detectors are more sensitive to the phase of the wave than its amplitude when one carries out param- 
eter estimation by matched-filtering the data with a template model. Because the covariance matrix of the noise is diagonal 
in the Fourier domain, one typically carries out matched-filtered parameter estimation in frequency space. The Fourier 
transform of the waveform for inspiraling neutron stars contains the term**-47 5W(f) œ [f(1)As +8(n)Aa] (amf), 
where f(n) and g(7) are functions of the symmetric mass ratio n = mim /m?, while Às a = (A1 + Az) /2 are the non- 
dimensional, symmetric and asymmetric combinations of the tidal deformabilities with Ay = A,/ mÀ representing the 
dimensionless tidal deformability for bodies A = (1,2). Extracting both combinations from the data from this single term 
in the Fourier phase seems impossible due to degeneracies among them, but there are two ways out??. 

One option is to choose an equation-of-state model to compute Asa as a function of the mass of the stars and determine 
the posteriors for the parameters of the model equation of state and the central densities by comparing to the data. Another 
option is to use equation-of-state insensitive relations. In the latter approach, one uses the binary Love relations*®*’ to 
prescribe A, in terms of A, analytically, thus making OW (f) a function of only one A; and n. One can then carry out 
parameter estimation to extract both A, and n (because the symmetric mass ratio also appears in other terms of the Fourier 
phase independent of the tidal deformabilities). Once A, has been extracted, one can then use the binary Love relations 
again to extract ae and from knowledge of both of these combinations one can trivially extract both A, and Az, without 
ever choosing an equation-of-state model. 


A caveat of this approach is that for the binary Love relation to apply, the two neutron stars in a binary must both be on the 
primary stable branch; for example, the relation is inapplicable for twin stars if one star is on one stable branch and the other is 
on a different, higher-density stable branch>». Care also needs to be applied when using the binary Love relation to place lower 
bounds on the tidal deformability, because these lower bounds can be below what is realistic for neutron stars and therefore 
may be outside the region in which the binary Love relations are valid>°. 

Both approaches have been applied by the LIGO/Virgo Collaborations on GW170817°* and the posteriors of the mass and 
radius are shown in Fig. 4 (see also Ref.*® for related work). Not surprisingly, both approaches lead to consistent posteriors in 
the mass-radius plane in the sense that the 50% credible regions overlap between the two analyses, suggesting that indeed the 
data is more informative than any systematic error incurred in the modeling. 


5 The power of coincidence 


But there is more! The double neutron star coalescence GW170817, which we recall occurred just 40 Mpc away, had 
counterparts over the entire electromagnetic spectrum>®. Gamma rays indicative of a gamma-ray burst observed ~ 20 — 30° 
off-axis were seen just ~ 1.7 seconds after the peak of the gravitational wave event. This was followed by ultraviolet through 
infrared emission over hours to weeks, with X-rays first detected nine days after the initial event and radio waves visible even 
now. 

The overall picture is consistent with predictions made prior to the event, in which the energy released by the coalescence 
of the neutron stars emerges in multiple forms: (1) gravitational waves, (2) a short gamma-ray burst (in which the gamma 
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Figure 4. Posterior for the mass and radius of each (blue and orange) components of the binary system that generated 

GW 170817, using a parametrized equation-of-state model (left panel) and the universal-relations method (right panel). The 
solid and dashed boundaries in each panel correspond to the 90% and the 50% credible region. A few mass-radius curves are 
overplotted in gray for comparison. The solid lines of the top left represent a non-rotating black hole, whose radius is twice its 
mass (times G/c? for unit consistency), and the so-called Buchdahl limit?” (ie. the maximum gravitational compactness that 
static and spherically-symmetric matter configurations must satisfy in general relativity under certain conditions on the energy 
and pressure). This figure is taken from Ref.°®. 


rays are produced by a blast wave moving with a Lorentz factor of at least hundreds), and (3) far more prolonged emission 
produced by the quasi-spherical, and much slower (v ~ 0.01 — 0.1c) outflow of unbound matter that was in the neutron stars 
(where the radiated energy comes from the decay of heavy radioactive nuclei). This last component, which might have had 
~ 0.01 —0.1 Mo in total mass, has been dubbed a “kilonova" or “macronova"*’. The highly neutron-rich outflow is thought to 
produce heavy elements, such as lanthanides and actinides, with high efficiency; thus, double neutron star mergers, and possibly 
mergers between black holes and neutron stars, could produce most of the heavy elements in the universe. 

GW 170817 and events like it provide us with ways to place an upper limit on the maximum mass Mmax of a non-rotating 
neutron star, albeit with astrophysical caveats®”°?. The argument is that if after the merger of two neutron stars the remnant is 
a long-lived and rapidly rotating neutron star, then if the process of merger generated a strong magnetic field the star would 
slow down rapidly by magnetic braking. This would then inject a large fraction of the ~ 10° erg rotational energy into the 
afterglow and kilonova. This excess energy is not observed in GW170817, which suggests that instead the merger remnant 
rapidly collapsed to a black hole. As a result, the total mass of the binary had to be larger than what can be supported by 
a rigidly rotating neutron star, and thus the maximum mass is bounded from above. The caveats are that there is no direct 
evidence for the formation of a black hole (as there would be if the gravitational wave data were many times more precise than 
they were), and the production of magnetic fields to slow the star’s rotation is highly uncertain. 

Nonetheless, under these assumptions, the estimated total mass of the double neutron star binary and the assumption that the 
remnant collapsed quickly imply that Mmax for a nonrotating neutron star is less than ~ 2.2 Mo; note that this relies on a fairly 
well-understood translation between the maximum mass of a rotating neutron star, such as is formed in the merger, and the 
maximum mass of a nonrotating star. Detailed numerical models of the outflow as inferred from electromagnetic observations 
yield similar answers, with implied values for Mmax in the range ~ 2.2 —2.3 Mor. 

If these upper limits are reliable, then, in concert with the existence of a few ~ 2 Mo neutron stars, they provide a remarkably 
tight constraint on Mmax: just ~ 2 — 2.3 Mọ to be conservative. This would eliminate both soft equations of state (which 
have Mmax < 2 Mo) and hard equations of state (which have Mmax > 2.3 Mo), therefore cutting down considerably on viable 
descriptions of the dense matter inside neutron stars. The most important addition to this information is independent, reliable 
and precise measurements of neutron star radii, which we now discuss. 
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6 A NICER way to nuclear astrophysics 


Precise neutron star radii have long been coveted by nuclear physicists, because they would arguably discriminate between 
different equation of state models better than any other single measurement. As a result, numerous radii have been reported over 
the years, usually focused on X-ray observations of spectra integrated over many neutron star rotation periods. It has, however, 
become evident that this method is susceptible to potentially serious systematic errors. For example, the X-ray emission from 
a typical cooling neutron star without pulsations can be fit equally well using a pure hydrogen atmosphere, a pure helium 
atmosphere, and even a black body spectrum (although neutron stars do not emit as black bodies). Despite the equally good fits, 
the inferred radii differ dramatically depending on the assumptions; for example, hydrogen and helium atmospheres can give 
radii that differ by as much as 50%. Thus, even if the formal statistical precision of the radius measurement is excellent, the 
reliability may not be. 

The Neutron star Interior Composition Explorer (NICER) adds a new dimension to the X-ray observations. In additional 
to its other observing tasks, NICER has been pointed at a small set of non-accreting pulsars for more than a million seconds 
each. These pulsars have X-ray emitting “hot spots" rotating with the star on the stellar surface that are believed to be produced 
by the impact of high Lorentz factor electrons and positrons on the surface, which are generated as part of the process that 
produces radio pulsar emission. NICER records the arrival time of each photon to better than 100 nanosecond accuracy, which 
is much shorter than the ~few millisecond rotational periods of these pulsars. The data can therefore be considered as spectra 
as a function of rotational phase, which is why it is sometimes referred to as time-resolved X-ray spectroscopy. Although more 
studies need to be performed, existing work suggests that when rotational phase information as well as spectra are obtained, 
then if the fit to the data is statistically good, the inferred radius and mass will not be significantly biased. For example, 
although the models of the shape and temperature distribution of hot spots cannot be perfectly correct, use of such models will 
either (1) produce a statistically poor fit, which then motivates the development of better models prior to radius inference, or 
(2) produce a statistically good fit, in which case the inferred radius and mass can provisionally be accepted to be reliable as 
well as precise. 

Like inference from LIGO data, inference of the radius from NICER data proceeds along standard Bayesian lines: given 
a model for the time-dependent spectra with parameters and associated priors (e.g., the mass, radius, spot shapes, locations, 
and temperatures, and the observer inclination angle), the NICER team determines the likelihood of the data given the model 
with specific parameter values, and iterates using a sampler until they obtain the posterior. Different parameters affect the 
phase-dependent spectra in ways that can be partially degenerate. For example, weak modulation could be produced by a very 
small spot (whose flux might be less than the background flux), or a spot that covers almost the entire star, or a spot nearly 
centered on the rotational pole, or an observer inclination nearly aligned with the rotational pole. However, with the hundreds 
of thousands of counts obtained in NICER observations, these degeneracies can be broken. 

At present, two independent groups within the NICER team have inferred the mass and radius of two pulsars: PSR JO300+0451 
and PSR J0740+6620. The first one is an isolated pulsar (which, being isolated, does not have an independently measured mass) 
that spins at a frequency of 205.53 Hz. For PSR JO030+0451 the teams reported a mass of ~ 1.4 Mo and 68.3% credible regions 
on the radius of 12.0 — 14.3 km (Miller et al. 2019°’; see the left panel of Figure 5) and 11.5 — 13.9 km (Riley et al. 2019°). 
The teams also found very similar hot spot locations and shapes. These radii are compatible with the 90% credibility radius 
upper bound of ~ 13.5 km found from GW170817, which contained two neutron stars which also had masses ~ 1.3 — 1.4 Mo. 

The second pulsar, PSR JO740+6620, spins at 346.53 Hz and is ~ 20x fainter in X-rays than PSR JO300+0451, but it is ina 
binary system, so one possesses independent measurements of the mass and observer inclination angle from radio observations. 
The Green Bank and CHIME radio telescopes have inferred a mass of M = 2.08 + 0.07 Mo and an observer inclination to the 
binary orbital axis of Oops ~ 87.5°. The inclusion of this independent information is crucial to obtain a good measurement of 
the neutron star radius. In addition, XMM-Newton data for this pulsar produced a more precise estimate of the pulsar X-ray 
flux than was possible using NICER alone, given the low flux of the pulsar and the comparatively high X-ray background in 
NICER observations. Combining the NICER, the XMM-Newton, and the radio data, the two teams reported 68.3% credible 
regions on the radius of 12.2 — 16.3 km (Miller et al. 202168; see the right panel of Figure 5) and 11.4 — 13.7 km (Riley et 
al. 20217°). The difference in the credible regions is primarily due to the use of different statistical samplers and different 
assumptions about the cross-calibration between NICER and XMM-Newton. 

The net result of the NICER measurements is that the radii of ~ 1.4 Mo and ~ 2.1 Mọ neutron stars are not very different 
from each other (indeed, they are consistent with being the same), with both being on the order of 12 — 14 km. This implies a 
relatively hard equation of state, with a maximum mass sufficiently above the ~ 2.1 Mo mass of PSR JO740+6620 that the 
radius has not yet bent toward smaller values (which is characteristic of the mass-radius curve near the maximum mass). The 
tidal deformability measurement for GW170817 eliminates the hardest equations of state, so together NICER and gravitational 
wave measurements have significantly narrowed the plausible list of candidates for the high-density equation of state. 
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Figure 5. Left panel: posterior probability density for the mass and radius of the isolated neutron star PSR J0030+0451 using 
only NICER X-ray data (original figure from®’). Right panel: posterior probability density for the mass and radius of the binary 
neutron star PSR J0740+6620, using NICER and XMM-Newton X-ray data as well as Green Bank and CHIME radio data 
(original figure from®*), These results, and the consistent results from the independent analyses off? and”?, imply that the 
radius of a neutron star is roughly constant from ~ 1.4 Mọ to ~ 2.1 Mo, and thus that matter in the core of neutron stars has 
relatively high pressure. 


7 The beauty of future dreams 


In the last 6 years, we have gone from living in a metaphorical data desert, with only a handful of gravitational wave observations, 
to living in a metaphorical data forest, with over 50 observations and counting. Clearly, gravitational wave astrophysics is here 
to stay, and we have only but began to explore this forest. The fourth observing (04) run of the LIGO and Virgo detectors, 
which is expected to reach design sensitivity, is scheduled to start in 2022, and the fifth observing (O5) run, which is expected 
to reach even better sensitivities, will start a few years later. The observing range for O4 will be about 50%—90% larger than 
that of O3, and the range for O5 is expected to be three times larger than that of O3. Since at the low redshifts relevant for 
double neutron star observations the accessible volume scales with the cube of the observing range, one can expect many, many 
more observations of binary black holes, binary neutron stars and mixed binaries, perhaps even in the hundreds by the time of 
O5 in just 1 year of data. By the mid-2030s, third generation ground-based detectors such as the Einstein Telescope and the 
Cosmic explorer will improve sensitivity by an additional factor of 10 — 20. 

What will the advent of such an increase in sensitivity do for us? As we expect to detect many more events, we in particular 
expect more binary neutron star observations. The low signal-to-noise ratio binary neutron star events may be too far away to 
lead to an electromagnetic counterpart, but if some of these mergers occur close to the Milky Way (say at tens of Mpc), then 
there is a good chance of aGW170817 repeat. This time, however, since the gravitational wave sensitivity will be significantly 
larger, the signal-to-noise ratio for an event at 40Mpc could be three times larger than that of GW170817, and thus, close to 
100. Such a high signal-to-noise ratio gravitational wave event would inaugurate the era of precision gravitational wave nuclear 
astrophysics, since it would allow for a much more accurate measurement of the tidal deformabilities, and thus of the radius of 
neutron stars. In fact, for such a loud event, not only could one extract nuclear physics information from the inspiral phase of 
the event, but one may also be able to extract information from the merger and post-merger itself, which may be useful to probe 
e.g. the presence of a quark matter core inside hybrid stars’!. 

Meanwhile, NICER observations will continue, and in particular, more and more data will be accumulated of the (soon to 
be) three pulsars that have already been observed. In fact, in 2022 the NICER team is expected to publish the measurement of 
the radius of the third pulsar, and to update the estimated radii of PSR JO030+0451 and PSR JO740+6620 using new data and a 
better characterization of the instrument. These measurements will improve the precision of our knowledge of the sizes of 
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neutron stars over a factor of 1.5 in mass, which will therefore yield additional valuable information about the equation of 
state of matter beyond nuclear density and the existence of quark matter”? inside neutron star cores. Beyond NICER there are 
numerous planned and proposed facilities and missions that will extend our reach greatly: these include radio observatories 
(the Square Kilometer Array and the next generation Very Large Array) and numerous X-ray missions (such as Athena, the 
enhanced X-ray Timing and Polarimetry mission [eXTP], and STROBE—X), which can probe the masses, radii, moments of 
inertia, and cooling properties of neutron stars. 

Whatever the future may hold, what is clear is that the combination of information from electromagnetic and gravitational 
wave observations is revealing the states of matter that are realized in the cores of neutron stars, at wonderfully large pressures 
and densities. In fact, it is not unreasonable to wager that within the next 10 years, the equation of state of matter at a few times 
nuclear saturation density will be, for the first time, constrained to better than 10% in the low-temperature neutron star region of 
the quantum chromodynamics phase space. The challenge now and in the near future will be to find ever more creative ways to 
connect these observations to fundamental nuclear physics. The future, is thus, beautifully exciting. 
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